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1  SUMMARY 


This  document  describes  the  work  performed  as  part  of  the  Ananke  simulation  development.  This 
work  includes  the  addition  of  higher-fidelity  models  in  CU-TurboProp  and  validation  of  the  new 
tools.  The  description  provided  here  includes  the  mathematical  background  and  description  of 
the  models  implemented,  as  well  as  a  definition  of  the  software  interface  for  the  larger  Ananke 
tool.  Comparisons  of  the  software  with  semi-analytic  techniques  and  other  available  high-fidelity 
propagation  tools  are  also  presented,  and  all  tools  met  all  required  and  expected  accuracies. 

2  INTRODUCTION 

This  document  describes  the  work  performed  by  the  Colorado  Center  for  Astrodynamics  Research 
(CCAR)  at  the  University  of  Colorado  Boulder  as  part  of  the  Space-Based  Search,  Detect,  Track 
project.  This  project  assembled  tools  developed  by  several  organizations  to  demonstrate  the  feasi¬ 
bility  of  resident  space  object  (RSO)  characterization  using  optical  observations.  This  simulation, 
called  Ananke,  provides  an  end-to-end  demonstration  of  the  simulation  and  processing  of  opti¬ 
cal  observations  of  a  previously  unknown  RSO  with  the  final  output  describing  likely  properties 
of  the  object.  For  example,  the  integrated  tool  may  be  used  to  identify  an  object  as  an  active, 
nadir-pointing  spacecraft,  and  includes  a  probability  of  this  result. 

As  part  of  the  Ananke  simulation,  CCAR  provided  the  high-fidelity  six  degree-of-freedom  (6- 
DOF)  models  used  for  simulation  of  the  data  and  in  the  various  filters  used  to  profile  the  object. 
This  document  provides  an  overview  of  these  tools  and  describes  the  work  performed  as  part  of 
the  Ananke  project.  The  next  section  provides  an  overview  of  CU-TurboProp  and  its  role  in  the 
Ananke  simulation.  Section  2.2  then  summarizes  the  work  performed  over  the  one-year  period 
of  performance.  Section  3  describes  the  models  added  as  part  of  Ananke,  while  Section  3.12 
provides  an  interface  description  for  the  CU-TurboProp  software  in  Anake.  Section  4  then  presents 
validation  results  of  software,  followed  by  a  summary  and  description  of  future  work. 

2.1  CU-TurboProp  for  Ananke 

First  developed  in  2006,  CU-TurboProp  (hereafter  referred  to  as  TurboProp)  is  the  principle  tool  for 
rapid  orbit  propagation  at  CCAR.  This  software  provides  a  library  of  functions  to  quickly  propagate 
precise  orbital  or  surface  trajectories  in  the  convenience  of  the  MATLAB  and  Python  programming 
environments.  To  speed  up  the  computations,  the  orbit  propagators  use  initial  value  ODE  solvers 
coded  in  the  C  programming  language.  The  orbit  propagators  may  be  called  as  MATLAB  functions 
through  a  MEX  interface  or  in  Python  via  the  Simplified  Wrapper  and  Interface  Generator  (SWIG) 
module.  Hence,  the  user  gets  both  the  superior  execution  speed  of  C  code  along  with  the  ease  of 
matrix  computations  and  plotting  that  come  with  high  level,  interpreted  languages. 

TurboProp  was  originally  written  to  propagate  orbits  in  the  Circular  Restricted  Three-body 
Problem  and  convert  them  to  the  Jet  Propulsion  Laboratory  (JPL)  ephemeris.  Since  then,  and  par¬ 
tially  as  a  part  of  the  Ananke  project,  capabilities  for  higher- fidelity,  Earth-centered  scenarios  have 
been  added.  This  has  enabled  space  situational  awareness  (SSA)  research  ongoing  at  CCAR.  This 
research  first  came  in  the  form  of  testing  the  prototype  cubed-sphere  model  [1,  2],  which  was  de¬ 
signed  for  space  surveillance  and  reduces  the  computation  requirements  of  the  gravity  perturbation 
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model.  Since  then,  and  partially  as  a  part  of  the  Ananke  simulation,  capabilities  for  Earth-centered 
orbits  have  been  expanded. 

For  the  Ananke  simulation,  CCAR  provided  a  reduced  set  of  TurboProp  functions  to:  (1) 
streamline  the  interface,  (2)  work  with  MATLAB  defined  integrators,  e.g.  ode 4  5  () ,  used  in 
Ananke,  and  (3)  provide  the  necessary  inputs  to  the  bidirectional  reflectance  distribution  function 
(BRDF)  model  provided  by  Pacific  Defense  Solutions  (PDS)  for  Ananke.  Figure  1  provides  a 
high-level  overview  of  this  interface. 


Ananke  Simulation 


Orbit  Propagator 


Figure  1.  Integration  of  TurboProp  with  Ananke 


2.2  Work  Performed 

As  part  of  Ananke,  CCAR  expanded  and  updated  the  suite  of  models  used  in  TurboProp.  Table  1 
illustrates  the  status  of  the  software  before  and  after  this  project.  The  improvements  to  the  software 
are  further  described  in  the  following  sections. 

3  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 

Orbit  propagation  incorporating  six-degree-of-freedom  (6-DOF)  modeling  requires  the  computa¬ 
tion  of  modeled  accelerations  and  torques.  In  this  new  version  of  TurboProp,  the  model  is 

r  =  +  f’GP+SET  +  f*3B  +  f’Drag  +  f’SRP  +  ^ERP  +  ^GR  (1) 

m  =  rn2B+GP+SET  +  rnSRp  (2) 

where  the  subscripts  indicate  the  type  of  acceleration  ( r )  or  torque  (m)  included  in  the  dynamics. 
The  models  implemented  in  TurboProp  include  the  central  body  term  in  the  form  of  the  two-body 
equation  (2B),  gravity  perturbations  with  a  solid  Earth  tide  correction  (GP+SET),  third-body  per¬ 
turbations  (3B),  atmospheric  drag  (Drag),  solar  radiation  pressure  (SRP),  Earth  radiation  pressure 
(ERP),  and  general  relativity  (GR).  Only  some  of  these  terms  are  included  in  the  torque  model. 
The  following  sections  provide  an  overview  of  the  models  used  and  developed  for  the  Ananke 
simulation.  Since  Ananke  only  considers  objects  near  geosynchronous  orbit,  we  omit  a  discussion 
of  atmospheric  drag. 
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Table  1.  Overview  of  Ananke-Derived  TurboProp  Improvements 


Capabilities  Before 

•  Translational  motion 

•  High-fidelity  static  gravity  field 

•  JPL  Design  Ephemeris 

•  Cannonball  Models  (drag  and  SRP) 

•  IAU76/FK6  GCRF/ITRF  Model 


Capabilities  After 

•  Shifted  to  IERS  2010  Conventions 

•  New/Updated  models 

-  6-DOF  (Attitude  Dynamics) 

-  Gravity  Gradient  Torques 

-  Satellite  Plate  Model 

-  SRP  Torques 

-  I AU 2006/IAU 2000 Aro6 
GCRF/ITRF 

-  Solid  Earth  Tides 

-  Earth  Albedo  Model 

-  General  Relativity  Correction 

•  Developed  more  extensive  software 
validation  process 

•  Improved  software  efficiency  for  MAT- 
FAB  interface 


Although  not  immediately  apparent  in  Equations  1  and  2,  translation  and  attitude  motion  is 
coupled  through  the  SRP  and  gravity  field.  The  SRP  acceleration  is  dependent  on  the  orienta¬ 
tion  of  the  satellite  relative  to  the  satellite-Sun  vector.  Hence,  it  depends  on  the  current  attitude 
state.  Additionally,  the  current  position  of  the  satellite  influences  the  gravity  perturbation  and  SRP 
models,  creating  a  dependence  on  the  translation  state  in  attitude  dynamics.  This  illustrates  how 
high-fidelity,  6-DOF  modeling  of  the  satellite  motion  is  required  for  highly  accurate  state  propa¬ 
gation. 

3.1  IERS  2010  Conventions 

As  indicated  in  Table  1,  TurboProp  was  updated  to  more  closely  adhere  to  the  International  Earth 
Rotation  and  Reference  System  (IERS)  2010  conventions  for  orbit  propagation  [3].  This  mainly  re¬ 
quired  an  update  of  the  transformation  between  the  Geocentric  Celestial  Reference  Frame  (GCRF) 
and  the  International  Terrestrial  Reference  Frame  (ITRF).  More  information  on  this  update  may 
be  found  in  Section  3.11.  Otherwise,  new  models  included  in  TurboProp  are  based  on  the  IERS 
conventions  (where  appropriate). 
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3.2  Central  Body  Forces 

The  two-body  equation  describing  the  the  gravitational  forces  due  to  the  central  body  is 

r2B  =  ~^r  (3) 

where  r  is  the  position  vector  of  the  spacecraft  with  respect  to  the  planet’s  center  of  mass,  r  is  the 
magnitude  of  that  vector,  and  /i  is  the  gravitational  parameter  of  the  central  body.  In  this  document 

r  —  [x  y  z\T  (4) 

is  a  position  vector  with  components  x,  y,  and  z. 


3.3  Central  Body  Gravity  Perturbations 

The  aspherical  gravity  field  is  modeled  using  U,  an  aspherical  potential  function  excluding  the 
two-body  term  [4], 


where 


n= 2  m= 0 


UiXiViZ)  —  ^  ^  ^  ^  ^  (  r  )  ^ n’m  ,y  I  {En,rnEm  T  Bn,mFm} 


X 


Em  —  EiEm-i  —  FiFm-i,  Eq  —  1,  E\  —  — 


Fm,  —  FiEm-i  +  EiFm-h  Fq  —  0,  F\  —  — . 


(5) 

(6a) 

(6b) 


We  note  this  form  is  defined  in  terms  of  the  Cartesian,  not  spherical,  coordinates.  Thus,  it  does 
not  contain  the  polar  singularity  present  in  the  classical  formulation  presented  in  [5,  6,  7],  and 
elsewhere.  R  is  the  radius  of  the  central  body,  //  is  the  gravitational  parameter  of  the  central  body, 
and  n  and  m  are  the  degree  and  order  of  the  spherical  harmonic  model,  respectively.  Cnrn  and 
Srhm  are  the  normalized  Stokes  coefficients  describing  the  magnitude  of  the  spherical  harmonic 
functions,  and  are  specific  to  the  gravity  field  model  being  used.  Anjn  \a  —  f]  is  a  normalized, 
derived  Legendre  function  computed  using  the  following  recursion: 


A),o[o/I 
A-ifi  [a] 

^1,1  [cfl 
Antn  [cc] 

A n,m  \B\ 


1 

at/3 

V3 


2  n  +  1 


2  n 

C%Bn,m-An—  \tm.  [o^] 


An—  l,n—  1  [oi]  Tl  1 

I Jr 


Bn—  1  ,r 


~An- 


n—  2,m  Q: 


m  <  n 


(7) 


where 


Bn.m. 


(2n  +  l)(2n-  1) 
(n  +  m)  (n  —  m) 


(8) 
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The  Dn  m  and  R  n,m  coefficients  are  not  a  function  of  a,  and  are  precomputed  at  software  initial- 

’  %-l,m 

ization  for  future  use.  This  formulation  of  the  derived  Legendre  functions  is  adapted  from  personal 
correspondence  with  Dr.  Bob  Schutz  of  the  University  of  Texas  Center  for  Space  Research  and 
from  his  book  Satellite  Dynamics:  Theory  and  Application  being  written  with  G.  Giacaglia.  This 
adaptation  to  the  formulation  in  Schutz’s  book  was  based  on  [8],  and  is  consistent  with  [4]. 

The  acceleration  in  body-fixed  coordinates,  r,  is  found  by  taking  the  gradient  of  the  potential 
U  in  Cartesian  coordinates.  As  derived  in  [4],  the  following  equations  summarize  the  method: 


where 


and 


OO  /  ^  \  Th  H 


CL  i  — 


Cl  2  — 


«3  — 


n=2  '  '  m—  1 

OO  /  ^  \  TL  H 

E  7  E  TflA  n^m  i  Cn^mF1 m—  l) 

n=2  '  '  m—  1 

OO  /  ^  \  Tl  H 

^  ^  f  j  ^  ^  n,m  (^7i,m-®m  “1“ 


71=2 

oo 


<24  = 


71=2 


R 


771=0 
71  n 


E  7  E  (n  +  m  +  l)An,m  (C n  mEm  T  SnmFm ) 


m=0 


^GP  = 


/  (Z  \  r 

CL\ 

\ 

(  ~a3  +  O 4  ) - 

&2 

l  vr  Jr 

_Oj3_ 

) 

r«,tn.  nn,m+i/nn,m 


(n  +  m  +  1  )(n  -  m)(2  -  <50m ) 


n 


n,m 


(n  +  m)! 

m)!(2  —  Som)(2n  +  1)' 


(9a) 

(9b) 

(9c) 

(9d) 

(9e) 

(10) 

(ID 


In  the  software,  T„  m  is  precomputed  and  stored  at  initialization. 

The  equations  of  motion  are  integrated  in  the  inertial  frame,  so  the  spacecraft  state  vector  must 
be  converted  from  the  inertial  frame  to  the  body-fixed  frame  before  computing  the  accelerations. 
After  computing  the  body-fixed  accelerations,  they  must  be  converted  back  into  the  inertial  frame. 
The  GCRF/ITRF  conversion  from  an  inertial  to  a  fixed  frame  is  described  in  detail  in  Section  3.11. 

When  propagating  state  deviations  in  the  classic  linearized  filters,  the  state  vector  is  augmented 
with  a  state  transition  matrix  (STM).  In  such  a  formulation,  the  STM  $  maps  a  given  deviation 
vector  Xq  at  time  t0  to  a  given  time  tk 


xk  =  $(tk,t0)x  0.  (12) 

The  state  transition  matrix  is  typically  propagated  simultaneously  with  the  state  using 

®(tk,t0)  =  A(tk)$(tk,t0)  (13) 

where 


A(t)  = 


dF(X,t) 


dX(t)  ’ 
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X  is  the  state  vector,  and  A  is  of  the  form 


dx 

dx 

dx 

dx 

dx 

dy 

dz 

dz 

dy 

dy 

dy 

dx 

dy 

dz 

dz 

dz 

dz 

(15) 

dx 

dy 

dz 

dz 

dz 

-  dx 

dz  - 

Although  Ananke  does  not  require  the  STM,  the  gravity  gradient  formulation  presented  in 
Section  3.9  will  require  A.  The  remainder  of  this  section  presents  the  mathematical  formulation 
of  this  matrix. 

To  include  the  two-body  acceleration  in  the  variational  equations  used  to  create  matrix  A,  the 
partials  of  the  two-body  portion  of  r  with  respect  to  r  are  computed  in  this  manner: 


x 2  1  xy  xz 


A-2B  ~ 


dr 

dr 


3y 

^3 


xy  y 2  1  yz 


(16) 


xz  yz  z2  1 

7"2  fp  2  j.2  g 

To  include  the  gravity  perturbations  in  the  matrix  A,  the  partial  of  the  gravity  field  portion  of  r 
with  respect  to  r  is  computed.  These  variational  equations  were  derived  in  flU,  and  are  provided 
here  as  a  reference.  The  resulting  matrix  is 


where 


^4-gp  — 


dr, 


GP 


dr 

y 


=  5  t-1  *] 


F  G 

~rT 

O 

1 

~rT~ 

G  M 

hT. 

+  [r  d] 

-1  0 

dT 

+ 


N  —  A 

-n 

Q 


-n 

—(N  +  A) 
R 


Q 

R 

-A 


z 

e  —  - 
r 

k  =  [0  0  1]T 
d  =  e[Q  R  0]T+[S  T  0]T 
A  =  (eo3  +  04) 

F  —  L  +  e(Me  +  2(T>  +  U3))  +  A 
G  =  -( Me  +  P  +  a3 ). 


(17) 


(18) 


(19) 

(20) 

(21) 

(22) 

(23) 

(24) 
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Additionally, 


L  = 

OO 

-  E 

71=2 

(RA 

Kr  J 

n  n 

)  Es  (n  +  m  +  1)  (n  +  m  +  2  )AntTn(Cr 

771=0 

7,777-^777  “1“  Sn^rnFrri^ 

(25) 

M  = 

OO 

-  E 

71=2 

(RA 

\r  j 

n  n 

J  ^  ^  -4t7,,777+2^  77,777  (C77, 777 Em  “1“  Sn^mF1 m/ 

777  =0 

i 

(26) 

N  = 

OO 

-  E 

71=2 

(RA 

\r  J 

77  77 

|  ^  ^  771(771  Anrn(CnrnEm_2  F  *St7,' 

777  =2 

777  Ejn — 2  ) 

(27) 

D  = 

OO 

-  E 

71=2 

(RA 

Kr  J 

77  77 

)  ^  ^  rn(m  An^rn{C n^mF m—2  *5n,i 

777  =2 

nErfYl  —  2) 

(28) 

P  = 

OO 

-  E 

71=2 

(RA 

Kr  J 

77  n 

)  ^  ^  Tfi l)z4njm_|_irnjm(C<rl5m£1r7i 

777  =0 

i  H“  ^77,777-fm) 

(29) 

Q  - 

OO 

-  E 

71=2 

(RA 

V  r  j 

77  n 

)  ^  ^  ^^AnrnJr\Y^nrn  (^71, 777 -S' 777—1  Snrn 

777=1 

.-^777-1) 

(30) 

R  = 

OO 

-  E 

71=2 

(RA 

\r  J 

77  77 

)  777-^4 77, 772-(-l  rn5m(*S77,?777,-E1m_l  C^77,777 

777=1 

,-^777-1) 

(31) 

S  = 

OO 

-  E 

71=2 

(RA 

\r  J 

77  77 

|  ^  ^ ~ t-  777/  — \^7TlAnm(CnmEm_i  + 

777=1 

*5n,777-^777— l) 

(32) 

T  = 

OO 

-  E 

77=2 

( R 3 

Kr ) 

77  n 

|  ^  ^  H-  771  +  1 )  777- z4n?m  -Em_i 

777=1 

^77,777-^777— 1 )  • 

(33) 

Note  the  equation  for  S  corrects  a  mistake  in  [4].  Finally, 


^ n,m  nn!m+2/nn,m 


( n  +  m+  1  )(n  +  m  +  2  )(n  —  m)(n  —  m  —  1)(2  —  Smt0) 


(34) 


Once  the  matrix  AGP  is  computed,  it  must  be  rotated  from  the  Earth-fixed  frame  to  the  inertial 
frame  before  being  added  to  A2b-  This  is  done  with  the  Earth-fixed  to  inertial  rotation  matrix  Tf 
described  in  Section  3.11.  Hence, 


T  —  A2 b  +  [Acp]j  —  A2b  +  [T^]T  [AGp]e  [T 'i  ■ 


(35) 


The  gravity  model  used  in  TurboProp  still  requires  a  definition  of  the  Stokes  coefficients  Cri/rn 
and  Sn>m.  These  coefficients  are  provided  through  a  text  file,  with  several  options  provided  with 
the  software.  For  Earth-centered  orbits,  we  provide  the  GGM03C  (GGM0  3C  .  sha)  [9],  EGM2008 
(tide-free:  EGM2 0 0 8t f  .  sha,  zero-tide:  EGM200  8zt .  sha)  [10],  EGM98  (EGM9  8  .  sha)  [11], 
JGM-3  (JGM3  .  sha)  [12],  and  WGS-84  (WGS8  4  .  sha)  [13]  gravity  models.  We  note  that  the 
values  of  //  and  R  used  in  the  spherical  harmonic  model  must  match  those  dictated  by  the  model. 

The  gravity  field  format  was  based  on  the  file  format  for  the  lunar  gravity  field  jgllOOkl.sha 
found  at 

http : / / pds -geo sciences . wust 1 . edu /miss ions /lunarp/ shadr . html 
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Table  2.  Header  Line  Format  for  Gravity  Field  Files 


C  format  string:  "  %le,  %le,  %le,  %d,  %d,  %d,  %le,  %le" 


Parameter  Description 

Format 

Units 

Reference  radius  of  the  planet 

Scientific  notation 

km 

Gravitational  parameter 

Scientific  notation 

km3/sec2 

Gravitational  parameter  uncertainty 

Scientific  notation 

km3/sec2  (discarded) 

Degree  of  model  field 

Integer 

N/A  (discarded) 

Order  of  model  field 

Integer 

N/A  (discarded) 

Normalization 

Integer 

0  =  unnormalized 

1  =  normalized 

Reference  longitude,  normally  0 

Scientific  notation 

degrees  (discarded) 

Reference  latitude,  normally  0 

Scientific  notation 

degrees  (discarded) 

The  first  line  is  a  header  and  contains  eight  parameters,  three  of  which  are  used  and  the  rest 
discarded.  The  parameters  are  shown  in  Table  2. 

Note  that  these  parameters  are  comma-delimited.  If  extras  .mu  was  not  input  by  the  user, 
the  gravitational  parameter  and  reference  radius  from  this  header  line  will  be  used  for  all  gravity 
calculations  for  the  central  body.  The  user  can  provide  extras  .mu  and  extras  .  radius  to 
override  the  parameters  in  the  gravity  header  line.  Since  these  values  are  dictated  by  the  model 
selected,  e.g.,  GGM02C,  this  input  should  only  be  used  to  change  the  length  or  time  units  to  match 
those  of  the  propagated  state. 

After  the  header  line,  “data  lines”  contain  the  gravity  field  coefficients.  These  coefficients  are 
also  comma-delimited  also,  and  are  described  in  Table  3. 

Table  3.  Data  Line  Format  for  Gravity  Field  Files 


C  format  string:  "  %d,  %d,  ‘ 

E>le,  %le,  %le,  %le" 

Parameter  Description 

Format 

Degree  index  n 

Integer 

Order  index  m 

Integer 

Coefficient  Cnjn 

Scientific  Notation 

Coefficient  SUjm 

Scientific  Notation 

Uncertainty  in  coefficient  Cn  m 

Scientific  Notation 

Uncertainty  in  coefficient  Sn>m 

Scientific  Notation 

3.4  Solid  Earth  Tides 

The  solid  Earth  tide  model  now  implemented  in  TurboProp  uses  the  formulation  presented  in  the 
IERS  2010  conventions  [3].  The  SET  perturbations  account  for  temporal  variations  in  the  mass 
distribution  of  the  solid  Earth  due  to  other  planetary  bodies  (Moon,  Sun,  etc.).  It  does  not  account 
for  variations  in  the  liquid  distribution  (oceans,  etc.),  which  would  be  handled  by  a  liquid  Earth  tide 
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model.  The  theory  behind  such  models  is  extensive,  and  beyond  the  scope  of  this  document.  In  the 
interest  of  brevity,  we  only  provide  a  high-level  overview  of  the  SET  model  with  key  equations  to 
aid  in  understanding  the  elements  of  its  computation. 

When  augmenting  the  spherical  harmonic  model  presented  in  Section  3.3  to  account  for  tem¬ 
poral  variations,  the  SET  model  simply  perturbs  the  Stokes  coefficients 


r  —  r  4-  Af 


n,m 5 


Q  —  9  4-  A  ^ 

un.m  un.m  '  ‘—xur 


n,ra? 


(36) 


where  the  deviations  ACn  m  and  ASn  rn  are  determined  by  the  model.  The  acceleration  vector 
v gp+set  and  the  variational  equations  A2B+gp+set  are  found  by  evaluating  Equations  9e  and  1 8  with 
the  corrected  coefficients.  The  temporal  changes  are  modeled  using  several  terms,  including  varia¬ 
tions  in  the  pole  of  rotation,  secular  variations  in  the  low-degree  coefficients,  frequency-dependent 
effects,  and  others.  Secular  variations,  although  not  strictly  a  part  of  the  SET  model,  are  modeled 
using 

Cnfl(t)  =  Cnfl(to)  +  Cnfi(t  —  to)  (37) 

where  the  secular  rates  C'n,o  and  epoch  times  f0  are  provided  in  Table  4  (adapted  from  Table  6.2  of 
[3]).  The  remainder  of  the  current  SET  model  is  a  three  step  process: 

1 .  Evaluate  the  frequency-independent  terms  modifying  the  (n,  m)  coefficients  for  n  =  0, . . . ,  4. 

2.  Generate  corrections  due  to  the  frequency-dependent  contributions  in  the  (2,  m),  m  =  0, 1,  2 
terms. 

3.  If  required,  remove  the  permanent  tide  correction. 

After  computing  the  corrections  A Cn'jm  and  A  S$m  for  steps  i  =  1,2,3,  the  final  correction  is 
simply  the  sum  over  i. 


Table  4.  Secular  Gravity  Field  Variations  [3] 


Coefficient 

Value  at  J2000.0 

Rate  /  yr  1 

6^2,0 

-0.48416948  xlO”3 

11.6xl0"12 

Qi;o 

0.9571612xl0~6 

4.9xl0"12 

<4,0 

0.5399659  xl0“6 

4.7  x  10“12 

The  Step  1  corrections  depend  on  the  position  of  any  perturbing  third  bodies,  usually  the  Sun 
and  Moon.  Mathematically,  the  corrections  A Cn,L  and  ASnX  for  the  n  —  2,3  terms  are 


A  Ck% 
A  s£L 


^12,171 

2  n  +  1 


Pn,m  Sin  (f>p 


cos(mAp) 

sin(mAp) 


(38) 


where  kn^m  is  the  nominal  Love  number  for  degree  n  and  order  m,  rp  is  the  distance  from  the  Earth 
to  the  third  body  (computed  using  the  JPL  design  ephemeris),  nP  is  the  third  body  gravitation 
parameter,  Pn.rn  is  the  associated  Legendre  function,  op  and  \p  are  the  geocentric  latitude  and 
longitude  of  the  third  body,  and  p  is  the  index  of  summation  over  the  third-bodies  considered  in 
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the  model.  The  ‘©’  subscript  indicates  a  parameter  for  the  Earth.  For  the  n 
coefficients, 


AG*1) 

4  ,m 

A  9O) 

^^4,m 


cos(mAp) 

sin(mAp) 


4,  m 


0,1,2 

(39) 


where  accounts  for  anelasticity  of  the  Earth.  The  values  for  kn  rn  and  k^n  may  be  found  in 
Table  6.3  of  [3], 

In  Step  2  of  the  SET  model,  we  generate  the  frequency-dependent  contributions.  For  the  (2,  0) 
term 

AC'S  =  (-^/(2.0)  cos(6)/(2,o))  -  0/(2, o)  sin(6l/(2,0)))  (40) 

/( 2,0) 


where  //(2,o)  and  Oja, a,  are  the  in-phase  and  out-of-phase  amplitudes  for  each  of  the  tides  /( 2,  0). 
Example  /( 2,  0)  tides  include,  but  are  not  limited  to,  the  Sa  and  Mm  tides.  0/(2, o)  is  defined  on 
p.  84  of  [3],  which  is  a  function  of  the  Greenwich  Mean  Sidereal  time,  the  fundamental  arguments 
of  nutation  theory,  and  the  multipliers  of  these  arguments  for  the  tide  /( 2,  0).  A  list  of  the  tides, 
the  amplitudes,  and  the  vector  multipliers  for  the  (2, 0)  correction  may  be  found  in  Table  6.5b  of 
[3].  Using  the  values  for  the  amplitudes  given  in  Tables  6.5a  and  6.5c  of  [3],  the  corrections  to  the 
(2, 1)  and  (2,  2)  terms  are 


1  1 

o?  rl 

<  < 

1 _ 1 

=  E 

T (2,1)  L 

<  < 

1 _ 1 

=  E  ' 

/( 2,2)  L 

2,1)  Sin(0/(2,l))  +  0/(2, 1)  COS(0/(2,1) 
^/(2,1)  cos(0/(2,i))  -  0/(2, 1)  sin(0/(2ji)) 

7/(2, 2)  COs(0/(2,2)) 

T  •_  //I 


(41) 

(42) 


Finally,  Step  3  removes  the  permanent  tide  correction  from  the  O2,o  term,  if  required.  Zero 
tide  models  of  the  gravity  field  include  a  correction  for  the  permanent  tide  already  accounted  for 
in  Steps  (1)  and  (2).  Hence,  this  correction  must  be  removed  to  prevent  its  duplication.  This 
correction  is  dependent  on  the  model,  and  is  not  necessary  when  using  a  tide  free  model. 

Although  not  described  here,  the  TurboProp  implementation  also  includes  corrections  for  the 
solid  and  ocean  pole  tides.  More  information  on  these  corrections  may  be  found  in  Sections  6.4 
and  6.5  of  [3]. 


3.5  JPL  Design  Ephemeris 

The  JPL  planetary  ephemeris  model  (DE)  expresses  the  location  of  third  bodies,  e.g.,  the  Moon  or 
Sun,  in  the  ICRF  inertial  reference  frame  centered  on  a  user-specified  central  body.  Upon  centering 
the  location  relative  to  the  Earth,  the  positions  are  expressed  in  the  GCRF  frame.  The  positions 
and  velocities  of  all  the  planets  are  obtained  from  JPL  DE  ephemerides:  DE403  [14],  DE405  [15], 
or  DE421  [16].  The  numeric  value  in  the  name  is  a  version  number,  with  DE421  being  the  most 
recent  publicly  available  version.  After  software  installation,  the  binary  ephemeris  files  for  DE403 
and  DE405  are  called  DE .  4  03  and  DE  .  4  04,  respectively.  Although  not  included  in  the  final 
release  of  TurboProp  for  the  Ananke  simulation,  CCAR  completed  the  integration  of  the  updated 
DE421  towards  the  end  of  the  Ananke  work  period. 
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The  code  used  to  create  the  ephemeris  files  and  to  perform  the  interpolation  was  adapted  from 
work  by  CCAR  graduate  student  Greg  Lehr  and  code  by  Mark  Hoffman  at 

ftp://ssd.jpl.nasa. gov /pub /eph /planets/ C-versions/hof fman/ 

The  ephemeris  data  files  were  obtained  from 


ftp ://ssd.jpl.nasa. gov /pub /eph /planets/ ascii 


Using  the  DE-determined  position  of  any  perturbing  bodies,  the  third-body  accelerations  are  then 
[5] 


^3  =  -/i3 


+ 


(43) 


where  //3  is  the  gravitational  parameter  of  the  third-body,  r:>  .s  is  the  vector  from  the  third-body  to 
the  satellite,  rc  3  is  the  vector  from  the  central  body  to  the  third-body,  and  r3  s  and  rc  3  are  the  vector 
magnitudes.  By  default,  the  gravitational  parameters  are  taken  from  the  JPL  ephemeris.  The  total 
perturbation  f3B  is  then  the  sum  of  the  accelerations  for  each  third-body  included  in  the  model. 


3.6  Earth  Radiation  Pressure  Model 


The  Earth  radiation  pressure  model  (ERP)  implemented  in  TurboProp  uses  the  formulation  pre¬ 
sented  in  [17].  Not  much  research  has  been  conducted  in  the  field  of  ERP,  and  this  model  appears 
to  be  the  one  most  commonly  employed  in  the  literature.  In  summary,  the  model  determines  the 
region  of  the  Earth  visible  to  the  satellite  as  a  function  of  its  Earth-centered  position.  The  visible 
area  is  then  segmented  into  N  elements,  each  with  some  surface  area  A,.  This  segmentation  is  il¬ 
lustrated  in  Figure  2.  The  acceleration  imparted  on  the  satellite  is  then  computed  for  each  element, 
with  the  sum  of  these  vectors  yielding  the  total  ERP  perturbation.  The  Earth  albedo  and  emis- 
sivity,  cij  and  eg  (for  element  j )  accounts  for  seasonal  variations  using  the  low-degree  Legendre 
polynomial  interpolation 


a,j  =  0.34  +  0.10  cos  (u(t  —  to))  Pi(sin(0j))  +  0.29P2(sin((/g))  (44) 

ej  =  0.68  —  0.07  cos(u(t  —  £0))Pi(sin((/g))  —  0.18P2(sin(</g))  (45) 

where  the  coefficients  are  the  empirical  values  presented  in  [17],  Pn (sin(/g ) )  is  the  n-th  degree 
Legendre  polynomial,  o3  is  the  equatorial  latitude  for  the  center  of  element  j,  uj  is  the  Earth’s  orbit 
mean  motion  («  27r/365.25  rad/yr),  and  t0  is  the  epoch  Dec  22,  1981.  With  these  parameters,  the 
ERP  acceleration  according  to  [17]  is  then 


where 


ri  = 


CrAs 

cMtc 


N 

U  ri’  =  X! 

3= 1 


Vo,/©  cos(9Sj)  +  ejE^ 


cos (og)Aj 


(46) 


(47) 


r  is  zero  or  one  if  the  center  of  the  element  is  in  shadow  or  illuminated,  respectively,  c  is  the  speed 
of  light,  IQ  is  the  local  solar  irradiance  for  the  element,  Cr,  M,  and  As  are  the  satellite’s  coefficient 
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of  reflectivity,  mass  and  area,  respectively,  Ee  is  the  exitance  of  the  Earth,  and  r  ?  is  the  element-to- 
satellite  distance.  We  note  that  the  presentation  here  expresses  the  equations  in  the  inertial  frame. 
For  visible  elements,  i.e.,  cos(ay  )  ^  0.0,  Equation  47  is  nonzero  if  the  element  is  not  illuminated. 
This  describes  the  force  imparted  on  the  satellite  due  to  continued  energy  emission  of  the  shaded 
portion  of  the  Earth. 


Figure  2.  Vector  geometry  of  the  TurboProp  ERP  model 


3.7  General  Relativity  Correction 

The  equations  for  gravitational  forces  discussed  in  previous  sections  obey  Newtonian  mechan¬ 
ics,  but  high-precision  propagation  requires  a  treatment  defined  by  the  theory  of  general  relativity 
(GR).  The  IERS  2010  conventions  define  a  correction  for  the  effects  of  GR  formulated  using  three 
principle  terms  [3,  Equation  10.12].  The  Schwarzschild  term  is  several  orders  of  magnitude  greater 
than  the  other  two.  Since  the  contribution  of  the  GR  correction  is  already  on  the  order  of  10-9  in 
low  Earth  orbit  when  compared  to  the  two-body  term,  we  truncate  the  smaller  terms.  This  is  further 
justified  since  the  relative  accuracy  of  the  gravity  field  does  not  exceed  10~9.  The  correction  for 
GR  in  TurboProp  is  then 


(48) 


3.8  Modeling  6-DOF  Motion 

High-fidelity  propagation  of  a  satellite  state  requires  properly  accounting  for  the  coupled  trans¬ 
lation  and  attitude  states.  To  this  end,  TurboProp  now  offers  a  capability  to  propagate  a  satellite 
state  with  these  additional  degrees  of  freedom.  TurboProp  uses  a  quaternion  (eft)  formulation  to 
represent  the  coordinate  rotation  from  the  inertial  ( i )  to  the  satellite  body  frame  (6),  and  includes 
propagation  of  the  angular  rates  based  on  a  given  torque  vector  (m).  This  section  describes  the 
formulation  of  the  attitude  dynamics. 
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Propagation  of  the  quaternion  is  modeled  as  a  composition  of  the  initial  quaternion  and  a  small 
rotation.  Given  the  initial  quaternion  q^, 

0  —Sou  i  —  Sou2  —S103 
5oj\  0  Sou^  —Su2 
Su2  —S0U3  0  ScUl 

(5cu3  S0J2  —Scui  0 

where  ou  is  the  pure  quaternion  generated  by  infinitesimal  angular  rate  vector  Sou.  We  note  that, 
for  notation  purposes,  we  include  the  scalar  portion  of  the  quaternion  first  as  qo,  followed  by  the 
vector  portion  q.  More  information  on  quaternions  may  be  found  in  [18]  or  [19]. 

Propagation  of  the  angular  rates  uses  the  Euler  rotational  equations  of  motion  [18] 

Jou  =  m  —  ou  x  J uo  (50) 

where  J  is  the  moment  of  inertia  tensor,  ou  is  the  rotation  vector  for  the  satellite  in  the  body  frame, 

and  m  is  the  vector  of  torques  in  the  body  frame.  The  vector  of  body  rotation  rates  ou  are  also 

expressed  in  the  body  frame.  Assuming  that  no  changes  to  the  satellite’s  mass  distribution  are 
modeled,  the  tensor  J*1  may  be  precomputed  at  initialization.  The  torques  computed  as  part  of 
TurboProp  are  described  in  the  following  sections. 

3.9  Gravity  Gradient  Torques 

The  previous  description  propagates  the  rotational  states  of  a  given  satellite  when  given  the  torques 
m.  In  TurboProp,  the  torques  can  include  those  due  to  gravity  gradient  effects.  The  full  central- 
body  gravity  model  uses  the  spherical  harmonic  model  previously  described  in  Section  3.3.  Also 
included  in  this  derivative  function  are  gravity  gradient  torques.  As  described  in  [4],  the  gravity 
gradient  torques  may  be  written  in  terms  of  the  variational  equations 

^■23(^33  —  J22)  —  A13J12  +  A12J13  —  ^23(^33  —  ^22) 

^2B+GP+SET  =  ^13(^11  —  J33)  +  A23J12  —  *4.12423  —  4i3(Ali  —  ^.33)  (51) 

_*4l2  (4^2  —  4]  1 )  —  A23J13  +  A13J23  ~  4i2  (*422  —  *4ll)_ 

where 

A  =  [-4-2b+gp+set] fj  =  n  (T?  [AGp+sET]e  Tl  +  Ai b)  Tf,  (52) 

i.e.,  the  gravity  field  variational  equations  (including  the  two-body  term)  in  the  satellite  body-fixed 
frame.  The  variables  Jmn  and  Anrn  refer  to  the  element  in  the  n-th  row  and  m-th  column  of  the 
J  and  A  matrices,  respectively.  The  matrices  AGP  and  AGP+SHT  were  previously  defined  in  Sec¬ 
tion  3.3.  The  rotation  matrix  TJ'  may  be  computed  from  q\,  and  T[  is  given  by  the  orientation 
parameters  (see  Sections  3.11  and  3.12).  We  note  that  Equation  51  includes  corrections  to  Equa¬ 
tion  8-30  of  [4].  The  software  in  the  Appendix  to  [4]  is  correct  and  matches  Equation  51  above. 

The  derivation  of  Equation  5 1  presented  in  [4]  requires  the  truncation  of  a  binomial  expansion. 
This  truncation  introduces  some  modeling  inaccuracy,  which  induces  a  small  error  when  compar¬ 
ing  the  implementation  to  simple,  completely  accurate  cases.  All  common  formulations  of  the 
gravity  gradient  torque  expressed  in  terms  of  J  include  this  error. 
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3.10  Shape-Dependent  SRP  Modeling 


Unlike  the  cannonball  model  for  solar  radiation  pressure,  the  model  provided  by  TurboProp  for 
Ananke  uses  a  satellite  plate  model  to  generate  an  attitude  dependent  acceleration  and  torque.  As 
presented  in  [20],  the  SRP  acceleration  may  be  expressed  as 


^SRP 


&k,i 


Ai 


-p(r)J2^cos(^kS^ 


k= 1 


(1  -  pk)rQs  +  2  (  —  +  pk  cos  fa  )  nk,i 


(53) 

(54) 


where  P(r)  is  the  pressure  at  r,  Ak,  pk,  and  8k  are  the  area,  specular  reflection  coefficient,  and 
diffuse  coefficient,  respectively,  for  the  k-\h  plate,  nkii  is  the  normal  vector  to  the  A  -th  plate  in  the 
inertial  frame,  and  f©s  is  the  unit  vector  from  the  center  of  the  plate  to  the  Sun.  Given  this  force, 
the  torque  is  then 

n 

ms RP  =  - P(r )  ^  Ak  cos  fa  (pkib  x  fask4)  ,  (55) 

k= 1 

where  pkj,  is  the  location  of  the  plate’s  center  in  the  body  frame.  Although  it  is  not  done  in  the 
current  implementation  due  to  time  constraints,  the  ERP  model  presented  in  Section  3.6  may  also 
be  formulated  using  this  shape-dependent  model. 

The  satellite  structure  is  provided  to  the  derivative  function  through  an  input  body  file.  This 
file  contains  information  on  the  plates  used  to  represent  its  shape  and  the  moment  of  inertia  ( J)  for 
the  satellite.  The  file  is  split  into  three  sections:  (a)  the  header,  (b)  the  moment  of  inertia  tensor  J 
and  (c)  the  individual  plates  that  constitute  the  model.  The  two-line  header,  described  in  Table  5, 
provides  metadata  on  the  file  and  the  total  mass  of  the  satellite.  The  second  section,  which  is  nine 
lines,  contains  the  elements  of  J.  A  description  of  the  file  format  for  this  section  may  be  found  in 
Table  6.  Finally,  the  remaining  lines  describe  the  locations  and  orientations  of  the  individual  plates 
in  the  model.  The  number  of  plates  was  provided  in  the  header,  and  Table  7  describes  the  format  of 
these  entries.  An  example  file,  ggSat  .bod,  is  provided  in  the  TurboProp/Data/  directory. 


Table  5.  Header  Format  for  Satellite  Body  Model  Files 


C  Format  Strings 

%d 

%e 


Parameter  Description 


Number  of  Plates 
Satellite  Total  Mass  (kg) 


3.11  GCRF/ITRF  Transformation 

This  section  describes  the  elements  that  must  be  considered  when  modeling  the  coordinate  trans¬ 
formation  between  the  Earth-fixed  frame  and  the  Earth-centered  inertial  frame.  The  discussion 
follows  the  one  presented  in  [21],  which  compared  the  effects  of  modeling  truncation  on  the  accu¬ 
racy  of  orbit  propagation. 
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Table  6.  Moment  of  Inertial  Line  Format  for  Satellite  Body  Model  Files 


C  format  string:  "%d  %d  %e" 

Parameter  Description 

Format 

Index  n 

Integer 

Index  m 

Integer 

Moment  of  Inertia  Tensor  Term  Jnm 

Scientific  Notation 

Table  7.  Data  Line  Format  for  Satellite  Body  Model  Files 


C  format  string:  "%e,  %e,  %e,  %e,  %e,  %e,  %e,  %e,  %e" 


Parameter  Description 

Format 

Plate  Area  (m) 

Specular  Coefficient 

Diffusion  Coefficient 

Plate  A"  Position  (m)  in  Satellite  Body  Frame 

Plate  Y  Position  (m)  in  Satellite  Body  Frame 

Plate  Z  Position  (m)  in  Satellite  Body  Frame 

Plate  Unit  Normal  Vector  Af-Component  in  Satellite  Body  Frame 
Plate  Unit  Normal  Vector  V'-Componcnt  in  Satellite  Body  Frame 
Plate  Unit  Normal  Vector  Z- Component  in  Satellite  Body  Frame 

Scientific  Notation 
Scientific  Notation 
Scientific  Notation 
Scientific  Notation 
Scientific  Notation 
Scientific  Notation 
Scientific  Notation 
Scientific  Notation 
Scientific  Notation 

Satellite  operations  requiring  ever-increasing  position  accuracy,  such  as  Earth  science  missions, 
are  on  the  rise.  Additionally,  software  utilizing  special  perturbations  for  orbit  propagation  is  rapidly 
replacing  routines  using  general  perturbations.  These  numerical  models  and  precise  position  so¬ 
lutions  require  accurate  frame  transformations,  particularly  between  the  Geocentric  Celestial  Ref¬ 
erence  Frame  (GCRF)  and  the  International  Terrestrial  Reference  Frame  (ITRF)  (i.e.,  inertial  i 
and  fixed  e).  The  procedure  and  models  needed  to  transform  between  the  GCRF  and  ITRF  are 
maintained  by  the  International  Earth  Rotation  Service  (IERS)  [3].  Earth  orientation  parameters 
(EOP)  are  crucial  for  correctly  executing  these  frame  transformations,  affecting  parameters  such  as 
station  coordinates,  satellite  positions,  and  non-spherical  gravitational  acceleration  vectors.  EOPs 
consist  of  the  following  [3]: 

•  Pole  Coordinates  ( xp,yp ):  coordinates  of  the  Celestial  Intermediate  Pole  (CIP)  with  respect 
to  the  IERS  Reference  Pole  (IRP)  in  the  International  Terrestrial  Reference  System  (ITRS). 
The  IRP  is  the  location  of  the  agreed  upon  terrestrial  pole  while  the  CIP  is  the  actual  axis  of 
Earth  rotation. 

•  Celestial  Pole  Offsets  (dX,dY):  observed  corrections  to  the  conventional  celestial  pole.  The 
conventional  celestial  pole  position  is  defined  by  the  IAU  Precession  and  Nutation  models. 

•  UT1  Time  Difference  (Af/Tl):  the  offset  of  Universal  Time  (UT1)  from  Universal  Coordi¬ 
nated  Time  (UTC),  i.e.  Af/Tl  =  UT1  —  UTC. 
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•  Length  of  Day  ( LOD ):  time  difference  between  the  observed  duration  of  a  mean  solar  day 
and  86,400  SI  secs. 

•  Atomic  Time  Offset  (A AT):  time  difference  between  International  Atomic  Time  (TAI)  and 
Universal  Coordinated  Time  (UTC),  i.e.  A  AT  =  TAI  —  UTC. 

Although  the  Atomic  Time  offset  is  not  technically  an  EOP,  it  is  included  in  this  list  because 
this  value  is  commonly  incorporated  into  tabulated  EOP  data  files.  EOPs  are  published  by  the 
IERS  once  per  day,  effective  at  0h  UTC  of  that  day.  These  daily  EOP  values  must  be  interpo¬ 
lated  for  maximum  accuracy  at  other  times  during  the  day.  In  literature  to  date,  however,  recom¬ 
mended  interpolation  methods  to  achieve  specific  accuracies  are  not  mentioned.  Additionally,  the 
published  values  of  the  pole  coordinates,  UT1  time  difference,  and  length  of  day  are  “tide-free”, 
meaning  that  users  requiring  high-precision  need  to  apply  additional  corrections  for  diurnal  and 
semi-diurnal  ocean  tides  and  libration  after  the  parameter  has  been  interpolated. [3]  This  fact  is 
fairly  unpublicized  within  astrodynamics  literature  [5,  6]. 

As  mentioned,  EOPs  are  needed  to  accurately  transform  from  the  ITRF  to  the  GCRF  and  vise 
versa.  Effective  Jan  1,  2009,  the  IERS  recommends  the  use  of  the  new  International  Astronomical 
Union  (IAU)  models  for  nutation  and  precession  to  perform  this  GCRF/ITRF  transformation. [3, 
22,  23,  24,  25,  26,  27]  The  new  model,  as  a  whole,  is  called  the  IAU  2000AR06/2006  model.  This 
includes  the  IAU  2000AR()6  nutation  theory  and  the  IAU  2006  precession  theory.  Several  methods 
for  implementing  this  new  model  have  been  developed  over  the  past  few  years,  but  are  all  within 
the  uncertainty  of  the  IAU  model  itself  [26,  25,  5,  28].  The  CIO-based  (Celestial  Intermediate 
Origin)  series  method  of  the  IAU  2000AR()6/2006  model  is  utilized  to  perform  the  GCRF/ITRF 
transformations  for  this  software.  Equations  (56-60)  below  describe  the  general  process  of  this 
transformation  and  more  details  on  the  IAU  model  and  implementation  processes  can  be  found  in 
[3,  5,  28], 


[w] 

[R\ 

[PN] 


ROT3(-s')ROT2(xp)ROTl(yp) 
ROT3  {-9era) 


1  -  aX2  -aXY 
-aXY  1  -  aY2 
-X  -Y  1 


X 

Y 

a{X2  +  Y2) 


ROT3(s) 


rGCRF  =  [PiV]  [R]  [W]  r  ITRF 
T gcrf  =  [PAf]  [R]  |  [W ]  rITRF  +  u>0  x  [W ]  rITRF 


For  ease  of  notation  in  other  sections  of  this  document, 


T[  =  [PA/]  [R]  [W: 


(56) 

(57) 

(58) 

(59) 

(60) 


(61) 


Figure  3  below  identifies  which  EOP  contributes  to  each  of  the  components  needed  for  the 
frame  transformation.  To  be  in  full  compliance  with  the  models  of  the  IERS  Conventions  (2010), 
namely  the  IAU  2000AR06/2006  model,  the  pole  coordinates,  UT1  time  difference,  and  length  of 
day  must  first  be  interpolated  to  the  appropriate  time  and  then  modified  by  ocean  tide  corrections 
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Figure  3.  Depiction  of  Which  EOP  is  used  to  Generate  Each  Element  in  the  Transformation 
Process. 


(■t'pj  Up)  (■t'pj  2/p) IERS  T  (Ax,  ^//)ocean  tides  (62a) 

Af/Tl  =  A(7T1iers  +  Af/Tlocean  tides  (62b) 

LOD  =  Z/OZ^iEJtg  A LO Doce&n  tides  (62c) 

The  libration  corrections,  (Ax,  Ay) Oration,  are  not  shown  in  Equation  (62)  and  are  not  included 
in  this  software  for  two  reasons:  1)  they  are  an  order  of  magnitude  smaller  than  the  ocean  tide 
corrections  and  2)  the  model  has  not  been  verified  by  observations  as  of  yet.  These  libration 
corrections  consist  of  diurnal  and  semi-diurnal  nutations  due  to  external  torque,  namely  luni-solar 
torque. [3]  The  corrections,  (Ax,  Ay)ocean  tides,  include  diurnal  and  semi-diurnal  variations  caused 
by  ocean  tides  for  the  three  EOPs  in  Equation  (62)  [3]. 

The  subscript  “IERS”  indicates  the  daily  published  values  that  have  been  interpolated.  These 
parameters  are  published  from  several  sources,  including  the  IERS,  National  Geospatial-Intelligence 
Agency  (NGA),  and  the  United  States  Naval  Observatory  (USNO).  However,  CelesTrak  has  com¬ 
piled  EOPs  from  all  available  sources  to  create  a  single  file  containing  EOP  data  from  1962  up  to 
predictions  1  year  into  the  future.1  The  file  maintained  by  CelesTrak  is  used  by  Analytical  Graph¬ 
ics  Inc.  Satellite  Tool  Kit  (STK),  and  TurboProp  requires  a  periodic  update  of  this  file.  Tools  are 
provided  for  this,  and  are  executed  as  part  of  the  installation  software. 

The  process  of  interpolating  EOPs  and  adding  corrections  is  outlined  by  the  flowchart  shown 
in  Figure  4.  It  is  prudent  to  remember  that  the  first  step  is  to  interpolate  any  EOP  of  interest.  The 
computation  of  the  ocean  tide  corrections  is  relatively  simple  and  only  requires  an  input  of  Julian 
Centuries  of  Terrestrial  Time,  Ttt ■  The  IERS  Conventions  provide  the  algorithm  for  computing 
the  ocean  tide  corrections  using  a  table  of  amplitudes  and  arguments  for  7 1  tidal  components  and 
FORTRAN  code  is  available  online  from  the  IERS  EOP  Product  Center  [3].  Two  different  source 
code  files  named  INTERP  .  F2  and  ORTHO_EOP  .  F3  are  available  online  for  the  computation  of 
ocean  tide  corrections.  The  file  ORTHO_EOP  .  F  is  based  on  the  full  model  from  [29]  and  the  file 

1  http : / /www. celestrak . com/ SpaceData/ 

2 INTERP  .  F  available  here:  ftp  :  /  /hpiers  .  obspm .  fr/eop-pc /models/ interp  .  f 

3  ORTHO-EOP  .  F  available  here:  ftp  :  /  /  tai  .  bipm .  org/ iers/conv2010/ chapter 8/ ORTHO_EOP  .  F 
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Figure  4.  Flowchart  of  EOF  Interpolation  and  Ocean  Tide  Correction  Process. 

The  INTERP  .  F  routine,  mentioned  above,  also  incorporates  EOP  interpolation  into  its  pro¬ 
cedure.  The  default  is  a  4-point  window  Lagrangian  interpolation  (i.e.,  3rd-order  Lagrange).  Al¬ 
though  this  is  the  default  interpolation  method,  it  does  not  mean  it  is  the  most  accurate.  [6]  suggests 
that  using  a  quadratic  polynomial  interpolation  for  the  pole  coordinates  and  Af/Tl  will  yield  suffi¬ 
cient  accuracy,  and  [30]  suggests  that  linear  interpolation  is  a  common  standard.  Various  methods 
of  interpolation  were  compared  in  [21],  with  this  implementation  leveraging  this  work. 

As  shown  in  Figure  5  below,  corrections  are  only  applied  to  the  pole  coordinates,  UT 1  time 
offset,  and  length  of  day  parameters.  The  celestial  pole  offsets  are  observed  corrections  to  the 
precession-nutation  theory.  The  corrections  are  oscillatory  as  would  be  expected  from  ocean  tides 
and  are  quite  distinct.  The  time  window  in  Figure  5  has  been  reduced  to  one  day  to  clearly  show 
the  effect  of  the  diurnal  and  semi-diurnal  ocean  tide  corrections.  The  corrections  to  Af/Tl  have 
a  RMS  amplitude  of  about  25  //sec  and  reach  as  high  as  73  /r sec  in  amplitude  during  the  3-month 
period  of  Jan  1  to  Apr  1,  2011. 

It  should  be  noted  that  high-rate  EOP  estimation  (sub-daily)  is  another  field  of  research  that 
is  ongoing.  The  models  described  and  used  in  TurboProp  are  constantly  in  the  process  of  being 
verified  and  improved  by  research  groups  performing  EOP  estimation  from  Global  Positioning 
System  (GPS)  and  Very-Long  Baseline  Interferometry  (VLBI)  measurements. [31,  32,  30]  It  is 
important  to  be  mindful  of  the  accuracy  of  each  model  and  keep  track  of  progress  being  made  on 
these  fronts,  because  it  directly  impacts  precise  frame  transformations  [21]. 
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Figure  5.  Ocean  tide  corrections  for  applicable  EOPs  on  19  Mar  2011. 


3.12  TurboProp/Anake  Interface 

As  mentioned  previously,  TurboProp  provides  an  interface  that  may  be  used  in  a  separate  MAT- 
LAB  integrator  for  the  modeling  of  the  principle  forces  and  torques  acting  on  a  satellite.  The 
software  provides  the  models,  written  in  C,  and  a  MEX  interface  that  may  be  used  to  access  the 
functions  from  MATLAB.  Extensive  details  on  the  inputs  may  be  found  in  the  CU-TurboProp 
manual  supplied  with  the  software. 

An  example  execution  of  the  model  has  been  provided  in  the  sample  code  SDT Sample .  m. 
The  principle  MATLAB  call  is: 

dY  =  derivmnex  (  'SixDOFJTull'  ,  time,  Y,  extras  ) 

where  time  is  the  time  of  the  input  state  Y,  the  extras  structure  provides  parameters  needed 
for  modeling  the  forces  and  torques,  and  dY  is  an  array  populated  by  Y  =  F(t,  Y).  Population  of 
the  extras,  re  f  t  ime  parameter  may  require  evaluation  of  the  UTC2  JED  ( )  function  to  convert 
to  ephemeris  time.  A  full  description  may  be  found  in  the  user’s  manual,  and  example  uses  are 
provided  in  the  sample  code. 

Also,  this  software  requires  internal  initialization,  i.e.  the  loading  of  several  data  files  describing 
the  force  model  environment,  etc.  The  current  version  uses  a  newer  implementation  that  maintains 
the  state  of  the  software  with  each  call  of  the  derivjnex  ( )  function.  To  clear  this  persistent 
memory  or  reinitialize  the  software,  the  user  must  first  call: 

clear  derivmnex 

to  clear  the  current  memory. 
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As  part  of  Ananke,  CCAR  developed  a  derivative  function  compatible  with  the  bidirectional 
reflectance  distribution  function  (BRDF)  software  developed  by  Pacific  Defense  Solutions,  LLC. 
The  BRDF  software  is  written  in  MATLAB,  and  cannot  easily  be  called  from  TurboProp  with¬ 
out:  (1)  porting  the  software  to  C,  or  (2)  calling  MATLAB  code  through  an  engine  running  on  a 
separate  thread.  Hence,  we  decided  to  not  directly  interface  TurboProp  with  the  BRDF  model. 
Instead,  within  the  MATLAB -level  derivative  function,  the  TurboProp  software  will  be  called 
first,  followed  by  the  BRDF-specific  function.  For  this,  we  define  a  new  derivative  function, 
SixDOF_Full_BRDF,  that  should  be  called  when  running  with  this  model.  Below  is  a  description 
of  the  interface. 

SixDOF_Full_BRDF  is  a  6-DOF  propagator  for  both  the  attitude  and  translation  state  of  a 
satellite,  and  includes  perturbations  from  the  aspherical  gravity  field,  atmospheric  drag,  and  third- 
body  effects.  The  modeled  torques  include  those  induced  by  gravity-gradient.  It  is  assumed  that 
solar  radiation  pressure  and  other  forces/torques  that  are  dependent  on  satellite  surface  properties 
will  be  handled  by  the  BRDF.  The  state,  yO,  is  a  14  x  1  vector: 


yO 

(  1) 

x  position  (km) 

yO 

(  2) 

y  position  (km) 

yO 

(  3) 

z  position  (km) 

yO 

(  4) 

x  velocity  (km/sec) 

yO 

(  5) 

y  velocity  (km/sec) 

yO 

(  6) 

i  velocity  (km/sec) 

yO 

(  7) 

q0  scalar  part  of  quaternion 

yO 

(  8) 

qx  component  of  quaternion  vector  portion 

yO 

(  9) 

q2  component  of  quaternion  vector  portion 

yO 

(10) 

q3  component  of  quaternion  vector  portion 

yO 

(11) 

uji  angular  velocity  in  body  1  axis  (rad/sec) 

yO 

(12) 

u>2  angular  velocity  in  body  2  axis  (rad/sec) 

yO 

(13) 

tu3  angular  velocity  in  body  3  axis  (rad/sec) 

yO 

(14) 

Cd  coefficient  of  drag  (non-dimensional) 

Position  and  velocity  coordinates  are  specified  in  the  GCRF  frame  centered  at  the  central  body, 
while  the  quaternion  q  describes  the  transformation  from  the  GCRF  frame  to  the  satellite  struc¬ 
ture  frame.  The  angular  velocity  vector  u>  describes  the  rotation  rate  of  the  body  relative  to  the 
inertial  frame,  and  is  expressed  in  the  structure  frame.  The  input  time  is  in  seconds  past  the 

extras  .  reftime. 

The  origin  of  the  coordinate  frame  is  at  the  center  of  mass  of  the  primary  body,  extras  .  mu 
must  be  included  and  should  contain  /i  (in  km3/sec2),  the  gravitational  parameter  of  the  primary 
body.  The  radius  of  the  primary  body  (in  km)  must  be  included  in  extras  .  radius.  Since 
this  model  employs  the  spherical  harmonic  model,  the  degree  and  order  of  the  expansion  is  in¬ 
cluded  in  extras  .  degord.  The  gravity  gradient  torques  may  include  perturbations  from  non- 
spherical  gravity  terms.  If  this  is  desired,  then  the  degree  and  order  of  the  variational  equations 
extras  .  degordstm  must  be  nonzero.  The  gravity  model  itself  (Cn  m  and  Sn/rn  coefficients) 
is  provided  in  extras  .  gravfi  le.  To  model  atmospheric  drag,  then  the  appropriate  values  for 
the  exponential  density  model  are  provide  in  extras  .  atmos.  In  this  case,  the  reference  density 
is  in  kg/m3,  and  the  scale  height  and  reference  radius  are  in  km.  The  area-to-mass  ratio  used  in 
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the  atmospheric  drag  calculations  are  in  m2/kg,  and  is  provided  to  the  software  in  extras  .  A_m. 
Set  extras  .  A_m  =  0  to  not  include  atmospheric  drag.  The  reference  time,  in  Julian  Ephemeris 
Time  (JED),  is  provided  to  the  function  in  extras  .  reftime.  Utilities  are  provide  in  the  MAT- 
LAB  tools  (see  TurboProp  users  manual)  to  facilitate  the  conversion  to  JED.  Planets  to  include  for 
third-body  perturbations  are  provided  in  the  extras  .planets  array,  while  the  central  body  is 
defined  in  extras  .  center.  If  the  central  body  is  specified  in  extras  .planets,  it  will  not 
be  included  in  the  calculation  of  third-body  perturbations.  The  positions  of  the  celestial  bodies, 
relative  to  the  central  body,  are  computed  using  the  JPL  DE  specified  in  extras  .  ephemf  lie. 
Finally,  the  satellite  moment  of  inertial  tensor  is  provided  in  ext  r  as  .  MO  I  as  a  3x3  matrix  in  the 
body  frame.  The  elements  in  extras  .MO I  must  be  in  units  ofkg-m2/rad2. 

This  version  of  the  software  also  requires  the  definition  of  several  files  to  define  the  rotation 
from  the  ITRF  and  the  GCRF.  These  new  elements  of  the  extras  structure  are: 

extras.XYsFile  File  containing  the  location  of  the  rotation  axis  in  the  X-Y  plane 
extras  .  EOPFile  File  containing  the  parameters  required  for  EOP  evaluation, 
extras  .  LeapSecFile  Definition  of  the  leap  seconds  used  in  time  conversions. 


A  sample  definition  may  be  found  in  the  sample  execution  files  provided  in  TurboProp. 

The  outputs  of  this  function  are  different  from  other  TurboProp  derivative  functions.  The  out¬ 
puts  are  a  vector  of  length  18,  which  are: 


out 

(  1) 

x  velocity  (km/sec) 

out 

(  2) 

y  velocity  (km/sec) 

out 

(  3) 

i  velocity  (km/sec) 

out 

(  4) 

x  acceleration  (km/sec2) 

out 

(  5) 

y  acceleration  (km/sec2) 

out 

(  6) 

5  acceleration  (km/sec2) 

out 

(  7) 

%  rate  of  change  for  scalar  part  of  quaternion 

out 

(  8) 

qi  rate  of  change  for  component  of  quaternion  vector  portion 

out 

(  9) 

q2  rate  of  change  for  component  of  quaternion  vector  portion 

out 

(10) 

q:>  rate  of  change  for  component  of  quaternion  vector  portion 

out 

(11) 

torque  about  body  1  axis  (N-m) 

out 

(12) 

torque  about  body  2  axis  (N-m) 

out 

(13) 

torque  about  body  3  axis  (N-m) 

out 

(14) 

rate  of  change  for  C'o  coefficient  of  drag  (non-dimensional,  should  be  zero) 

out 

(15) 

xsQ  satellite-to-sun  vector  in  inertial  frame  (km) 

out 

(16) 

ySQ  satellite-to-sun  vector  in  inertial  frame  (km) 

out 

(17) 

zsQ  satellite-to-sun  vector  in  inertial  frame  (km) 

out 

(18) 

v  fraction  of  Sun  visible  from  the  satellite,  1  =  fully  illuminated,  0  =  in  shadow 
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We  provide  an  example  interface  with  the  PDS  BRDF  model  in  the  file  exampleBRDF .  m 
in  the  top-level  directory  of  the  TurboProp  package.  Assuming  you  have  access  to  the  PDS  soft¬ 
ware,  then  you  can  execute  this  sample  code  using  runExampleBRDF  .m.  In  this  case,  change 
the  filePath  variable  to  point  to  the  directory  with  the  BRDF  test  cases  and  make  sure  the 
Multiple-Model  Adaptive  Estimator  (MMAE)  software  (developed  by  State  University  of  New 
York,  Buffalo)  is  in  the  MATLAB  path. 

4  RESULTS  AND  DISCUSSION 

4.1  Developed  Unit  Test  Cases 

The  TurboProp  software  includes  a  set  of  unit  test  cases  (UTCs)  to  exercise  various  components  of 
the  core  functionality.  The  unit  tests  currently  verify  the  following:  ITRF/GCRF  transformation, 
gravity  gradient  torques,  attitude  dynamics  and  solar  radiation  pressure. 

The  ITRF/GCRF  transformation  unit  tests  verify  the  EOP  and  CIP/CIO  parameter  lookup  and 
interpolation  function  as  expected  and  also  verify  the  end-to-end  transformation  code  operates  cor¬ 
rectly  (transformation  from  ITRF  to  GCRF  and  back  to  ITRF).  Since  the  transformation  software 
was  originally  developed  in  MATFAB  and  then  translated  to  C,  the  unit  tests  compare  the  MAT- 
FAB  and  C  results.  Furthermore  the  unit  tests  were  used  to  compare  the  C  implementation  of  the 
transformation  to  STK  output. 

TurboProp  generates  gravity  gradient  torques  from  the  spherical  harmonic  gravity  field  and 
spacecraft  inertia  tensor  as  described  in  Section  3.9.  An  alternative  method  models  a  simple 
spacecraft  as  a  collection  of  point  masses  and  calculates  the  gravity  gradient  torque  as  the  sum 
of  the  torques  from  each  point  [33].  The  unit  test  driver  computes  the  gravity  gradient  torques  via 
both  methods  for  10,000  random  spacecraft  locations,  and  verifies  the  maximum  relative  error  is 
10-6  N-m.  This  difference  is  expected  given  the  truncation  error  discussed  in  Section  3.9. 

To  verify  the  implementation  of  Euler’s  equations,  we  compare  the  propagation  to  analytic 
solutions  defined  under  torque-free  motion.  For  example,  when  given  an  axially  symmetric  body 
with  Jn  =  J-22  T33  (other  elements  zero),  then 

cos(o;p/) 
u }(t)  =  sm(u}pt) 

0 

where 

Up  =  0 Y —  -  1^  cc3.  (64) 

We  then  compare  the  numerically  propagated  00  using  Euler’s  equations  and  TurboProp  to  the 
solution  generated  using  the  this  analytic  solution.  This  case  is  tested  using  10,000  random  initial 
conditions  to  provide  more  confidence  in  the  result.  Tests  demonstrate  a  relative  accuracy  of  at  least 
10-13,  which  matches  the  expected  result  after  floating  point  error  accumulation  with  numerical 
integration  with  a  sufficiently  long  time  span. 

Unfortunately,  validation  of  the  SRP  and  satellite  shape  model  software  requires  more  com¬ 
plexity.  The  UTCs  developed  for  the  SRP  software  consider  two  models:  (1)  a  satellite  with 
12,800  faces  to  create  an  approximate  sphere,  and  (2)  a  flat  plate  with  the  normal  vector  parallel  to 


—  sin(a;pf)  0 

cos (upt)  0  u>(t0)  (63) 

0  1 
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rsat, ©•  In  the  first  case,  the  UTC  demonstrates  that  the  model  approximates  the  cannonball  model 
to  an  accuracy  of  10  m.  In  the  second  case  matches  the  cannonball  model  to  a  relative  accuracy 
of  10~13.  Both  of  these  accuracies  are  computed  over  10,000  random  initial  conditions.  In  the  first 
test,  the  orientation  of  the  modeled  satellite  is  randomly  varied.  In  test  (2),  the  location  of  the  plate 
in  Earth  orbit  is  randomly  chosen. 

4.2  Comparison  to  MMAE  Propagator 

In  the  process  of  integrating  the  TurboProp  software  with  the  Ananke  simulation,  we  also  compared 
to  the  default  propagator  included  with  the  MMAE  filter.  The  MMAE  filter  included  the  BRDF 
software  previously  mentioned,  and  provided  an  additional  comparison  to  ensure  the  TurboProp- 
BRDF  interface  functions  properly.  Through  this  process,  CCAR  assisted  in  improving  the  force 
models  included  in  the  MMAE  filter. 

4.3  Comparisons  to  Satellite  Tool  Kit  (STK) 

Satellite  Tool  Kit  (STK)  was  used  to  provide  independent  validation  of  various  TurboProp  com¬ 
ponents  such  as:  translational  state  propagation,  rotational  state  propagation,  spherical  harmonic 
gravity  field,  third  body  gravity,  JPL  ephemeris,  ITRF/GCRF  transformation  and  gravity  gradient 
torque.  However  it  is  not  possible  to  configure  STK  to  model  all  of  the  forces  TurboProp  includes, 
nor  is  it  always  possible  to  force  STK  to  model  some  forces  in  the  same  manner  as  TurboProp.  The 
primary  limitations  of  the  STK  validation  are  that  STK  does  not  couple  rotation  and  translation, 
STK  can  only  include  two-body  gravity  gradient  effects  and  it  is  difficult  to  mimic  the  TurboProp 
SRP  and  drag  models  in  STK. 

The  comparison  scenario  was  based  on  an  object  in  near  geostationary  orbit  with  the  body  and 
inertial  coordinates  initially  aligned.  The  object’s  initial  rotational  velocity  in  all  three  axes  was 
non-zero  and  the  start  of  the  simulation  was  Julian  day  2455197.516765046.  The  scenario  was 
propagated  for  seven  days  using  a  Dormand-Prince  8(7)  [34]  integrator  in  STK  and  TurboProp. 

The  force  model  for  both  STK  and  TurboProp  included  a  degree  and  order  30  EGM2008  tide 
free  spherical  harmonic  gravity  field,  two-body  gravity  gravity  gradient,  and  third  body  gravity 
from  the  Sun,  Moon  and  Jupiter.  The  force  model  did  not  include  drag,  SRP  or  translation/rotation 
coupling.  The  simulations  used  the  JPL  DE405  ephemeris  and  the  inertia  tensor  was  fully  popu¬ 
lated  (no  zero  elements)  with  a  mix  of  positive  and  negative  products  of  inertia. 

Figure  6  shows  the  difference  between  STK  and  TurboProp  for  both  the  translational  and  rota¬ 
tional  states.  The  translational  states  are  shown  in  radial,  in-track  and  cross-track  (RIC)  coordinates 
and  the  angular  velocity  is  expressed  in  the  body  frame.  The  maximum  RSS  position  error  over 
seven  days  was  25  mm  (approximately  10_8%  relative  error)  and  the  maximum  angular  velocity 
RSS  error  was  1.8  x  10-8  deg/sec  (approximately  10~3%  relative  error).  The  disagreement  in  po¬ 
sition  is  largely  due  to  the  inclusion  of  third  body  gravity  (the  error  drops  by  an  order  of  magnitude 
when  ignoring  third  body  gravity).  The  angular  velocity  is  different  in  that  the  increase  in  error  is 
not  linear;  rather,  it  increases  by  an  order  of  magnitude  approximately  every  2-3  days  (the  maxi¬ 
mum  relative  error  is  10~5%  after  1  day).  The  angular  velocity  difference  also  seems  to  be  more 
heavily  tied  to  integration  step  size  and  likely  derives  from  round-off  and  truncation  error. 
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Figure  6.  STK  vs.  TurboProp  Radial,  In-track  and  Cross-track  Differences  (left)  and  Angu¬ 
lar  Velocity  Differences  (right) 

5  CONCLUSIONS 

This  document  described  the  work  performed  as  part  of  the  Ananke  simulation  development.  This 
work  included  the  addition  of  higher-fidelity  models  in  CU-TurboProp  and  validation  of  the  new 
models,  each  of  which  were  described  in  the  document.  The  software  interface  for  the  Ananke 
simulation  tool  was  described,  and  validation  results  were  presented.  All  results  met  required  and 
expected  accuracies  in  the  tests  considered. 

Future  work  could  take  several  forms.  As  briefly  mentioned  in  the  overview  of  CU-TurboProp, 
the  software  includes  several  numerical  integration  routines.  Some  work  could  be  done  to  develop 
improved  methods  to  reduce  the  computation  time  required  for  orbit  propagation.  One  example 
includes  the  BLC-IRK  integrator  in  development  at  CCAR  [35,  36].  Additionally,  CU-TurboProp 
may  be  further  improved  by  adding  a  liquid  tide  model,  which  would  improve  accuracy  for  low- 
Earth  orbit  propagation,  and  the  modeling  of  geomagnetic  torques. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


A 

4 

-**-71,171 

BRDF 

r  C 

71,171')  u 71,771 

CCAR 
CIO 
CIP 
CU 
DE 
dX,  dY 
EGM 
EOP 
ERP 
F 

GCRF 

GR 

GGM 

GP 

GPS 

Jf 

IAU 

IERS 

IRP 

ITRF 

Jij 

J 

JGM 

JPL 

k 

rxjn,m 

LOD 

M 

m 

rn 

MMAE 

n 

Of 

ODE 

P 

p 

1  n,m 

PDS 

PM 

q 


Area 

Derived  Legendre  function  of  degree  n  and  order  m 
Bidirectional  Reflectance  Distribution  Function 
Stokes  coefficients 

Colorado  Center  for  Astrodynamics  Research 

Celestial  Intermediate  Origin 

Celestial  Intermediate  Pole 

University  of  Colorado 

Design  Ephemeris 

Celestial  pole  offsets 

Earth  Gravity  Model 

Earth  Orientation  Parameter 

Earth  Radiation  Pressure 

Force  model 

Geocentric  Celestial  Reference  Frame 
General  Relativity 
GRACE  Gravity  Model 
Gravity  Perturbations 
Global  Positioning  System 
In-phase  amplitude  of  tide  / 

International  Astronautical  Union 

International  Earth  Reference  System 

IERS  Reference  Pole 

International  Terrestrial  Reference  Frame 

i,  i  element  of  J 

Moment  of  inertia  tensor 

Joint  Gravity  Model 

Jet  Propulsion  Laboratory 

Geopotential  love  number  for  degree/order  ( n,m ) 

Length  of  Day 

Mass 

Spherical  harmonics  order 
Torque  vector 

Multiple-Model  Adaptive  Estimator 
Spherical  harmonics  degree 
Out-of-Phase  amplitude  of  tide  / 

Ordinary  Differential  Equation 
Pressure 

Associated  Legendre  function  of  degree  n  and  order  m 

Pacific  Defense  Solutions 

Earth  precession  and  nutation  matrix 

Quaternion 
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q  Vector  portion  of  the  quaternion 
g0  Scalar  portion  of  the  quaternion 
R  Earth  pole  rotation  matrix 
r  Magnitude  of  position  vector 
r  Position  vector 
r  Velocity  vector 
r  Acceleration  vector 
RIC  Radial,  In-Track,  Cross-Track 
RMS  Root  Mean  Square 
RSO  Resident  space  object 
RSS  Root  Sum  Square 
SSA  Space  Situational  Awareness 
SET  Solid  Earth  Tides 
SRP  Solar  Radiation  Pressure 
SWIG  Simplified  Wrapper  and  Interface  Generator 
6-DOF  Six  degree  of  freedom 
STK  Satellite  Tool  Kit 
STM  State  Transition  Matrix 

Tj  Transformation  matrix  from  frame  i  to  frame  j 
t  Time 

TAI  International  Atomic  Time 
3B  Third-body 
2B  Two-body 
U  Geopotential 

USNO  United  States  Naval  Observatory 
UTC  Universal  Coordinated  Time 
UT1  Universal  Time 
VLBI  Very-Long  Baseline  Interferometry 
W  Bias  matrix 
WGS  World  Geodetic  System 
x,  y,  z  Cartesian  position  components 
xp,  yp  Pole  coordinates 
x  State  deviation  vector 
X  State  vector 
5k  diffuse  reflection  coefficient 
A  Longitude 
u  Angular  rate 
//  Gravitation  parameter 
Pk  specular  reflection  coefficient 
p  Satellite  facet  location 
$  State  transition  matrix  (symbol) 

<f>  Geocentric  latitude 

<jl>  Angular  momentum  vector 

u  Pure  quaternion 
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